
clear
do "...\First.do"


********************************************************************************
* This do file uses patients who died prior to clinic closure to test for potentiel
* pre trend in the mortality outcome
* The test is done in several steps. 

** First, I identify all patient in the relevant clinics dying within 5 years of clinic closure
** Second, I predict new PCP (post closure if they were alive) from the surviving patients in the clinics
** Third, test if the predict new PCP SES explains mortality patterns before clinic closure
** This file does it on the clinic level
********************************************************************************

use "$work\analysis_sample.dta"
keep if inrange(age,40,70)


* Keep pnr x year of people in the analysis sample
preserve 
keep pnr year 
save "$work\sample_pnr.dta", replace
restore

* Keep PCP id (ydernr) of closing providers
keep if timing==-1
keep ydernr
duplicates drop
save "$work\temp_pre_ydernr.dta", replace



********************************************************************************
*** 1. Identify all patient in the relevant clinics dying within 5 years of clinic closure
********************************************************************************


********************************************************************************
* Load full pop
********************************************************************************
clear
use "$work\full_sample.dta"

* Keep relevant sample
keep if inrange(age,30,70)


********************************************************************************
** Keep people who dies
********************************************************************************
g help=(death==1)
bys pnr: egen keep=max(help)
keep if keep==1
drop keep

********************************************************************************
** Identify PCP one year before death
********************************************************************************
sort pnr year
g yder_0=ydernr if death[_n+1]==1 & pnr==pnr[_n+1] 

rename ydernr ydernr_copy
rename yder_0 ydernr

gsort pnr -ydernr 
replace ydernr=ydernr[_n-1] if pnr==pnr[_n-1] & ydernr==""

drop timing 
cap drop min_year


********************************************************************************
* Find patients who were in the clinincs that are in the analysis sample
********************************************************************************
merge m:1 ydernr using "$work\temp_pre_ydernr.dta"
keep if _merge==3 
drop _merge

* Merge closing year
merge m:1 ydernr using "$work\gp_closures_year.dta"
keep if _merge==3
drop _merge 

bys pnr: egen _year_max=min(year_max)
replace year_max=_year_max
drop _year_max

********************************************************************************
* Generate timing untill time of closure for the dying patients in the closing clinics
* keep 5 periods prior to clinic closure
********************************************************************************
 
g timing2=(year-year_max)
tab timing2 death, m

keep if inrange(timing2,-5,0)
tab timing2 death

********************************************************************************
* Keep only people who are NOT in the main analysis sample
********************************************************************************
merge 1:1 pnr year using "$work\sample_pnr.dta"
drop if _merge==2
g analysis_sample=(_merge==3)
drop _merge

bys pnr: egen max=max(analysis_sample)
drop if max==1

tab timing2 death

********************************************************************************
* Save dataset
********************************************************************************

save "$work\temp_all.dta", replace

********************************************************************************
********************************************************************************

clear 
use "$work\temp_all.dta"

keep pnr year timing2 death ydernr low_ses mean_age mean_male mean_dk solo ku au sdu other N_doctors male age non_dk married
keep if inrange(timing2,-5,0)

g pre_sample=1


********************************************************************************
* Second, predict new PCP (post closure) for dying patients from the surviving patients in the clinics
********************************************************************************

********************************************************************************
** Append with peer patients who are alive - they are used for the prediction
append using "$work\analysis_sample.dta"
********************************************************************************

********************************************************************************
* Timing variable
********************************************************************************
replace timing2=timing2 // Timing for pre-mortality patients
replace timing2=timing if timing!=. // Timing for non-pre-mortality patients
tab timing2 death, row

cap drop t2
g t2=timing2+4
tab timing2 t2
tab timing2 death if ydernr!=""

********************************************************************************
* GP identifier
********************************************************************************
g gp_fe2=gp_fe
replace gp_fe2=ydernr if pre_sample==1
drop ydernr 
rename gp_fe2 ydernr

********************************************************************************
* Keep relevant population
********************************************************************************
keep if inrange(age,40,70)

********************************************************************************
** Keep a balanced pre-period for surviving patients
********************************************************************************
cap drop help N_pre_2
g help=1 if inrange(timing2,-4,0)
bys pnr: egen N_pre_2=total(help)
drop if N_pre_2!=5 & pre_sample!=1 


********************************************************************************
* Merge ses of the "new SES " on this file is generated in XXX 
********************************************************************************
cap drop _merge
merge m:1 ydernr using  "$work\yder_year_mean_ses.dta", update
tab year _merge
drop if _merge==2
drop if _merge==1
drop _merge 

********************************************************************************
* Create quartiles - the share of surviving patients who gets a new low SES PCP
********************************************************************************
g copy=mean_ses_patient
sum mean_ses_patient if timing2==-4 & pre_sample==1, d
replace mean_ses_patient=. if !inrange(mean_ses_patient,`r(p5)', `r(p95)') // Take out extremes

cap drop _ses  ses
xtile _ses=mean_ses_patient if timing2==-4, n(4)
bys pnr: egen ses=max(_ses)

g _age0=age if timing2==-4
bys pnr: egen age0=max(_age0)


* Create dummies - in a certain quartiles or equal to or a higher order quartile
forvalues v=1/4 {
	
cap drop SES_`v'
g SES_`v'=(ses>=`v') if ses!=.
replace SES_`v'=0 if (ses<`v') & ses!=. 

cap drop SES_b`v'
g SES_b`v'=(ses==`v') if ses!=. 
replace SES_b`v'=0 if (ses!=`v') & ses!=.

}




********************************************************************************
* 3. test if the predict new PCP SES explains mortality patterns before clinic closure
** Focus on the pre-period [-4;0]
** Cap the age at time -4 to 64, as this matches the patients in the surviving sample used in the analysis
********************************************************************************


cap postclose john
postfile john T timing ses beta upper lower using "$table\PreDeath", replace


forvalues v=1/4 {
areg death ib4.t2##SES_b`v' male i.age married non_dk if low_ses==0 & inrange(t2,0,4) & age0<64, cluster(ydernr) a(ydernr)
mat A=r(table)
mat list A

post john (`v') (0) (0) (A[1,9]) (A[6,9]) (A[5,9])
post john (`v') (1) (0) (A[1,11]) (A[6,11]) (A[5,11])
post john (`v') (2) (0) (A[1,13]) (A[6,13]) (A[5,13])
post john (`v') (3) (0) (A[1,15]) (A[6,15]) (A[5,15])
post john (`v') (4) (0) (0) (.) (.)


areg death ib4.t2##SES_b`v' male i.age married non_dk if low_ses==1 & inrange(t2,0,4) & age0<64, cluster(ydernr) a(ydernr)
mat A=r(table)
mat list A

post john (`v') (0) (1) (A[1,9]) (A[6,9]) (A[5,9])
post john (`v') (1) (1) (A[1,11]) (A[6,11]) (A[5,11])
post john (`v') (2) (1) (A[1,13]) (A[6,13]) (A[5,13])
post john (`v') (3) (1) (A[1,15]) (A[6,15]) (A[5,15])
post john (`v') (4) (1) (0) (.) (.)


areg death ib4.t2##SES_b`v'##low_ses male i.age married non_dk if inrange(t2,0,4) & age0<64, cluster(ydernr) a(ydernr)
mat A=r(table)
mat list A

post john (`v') (0) (3) (A[1,37]) (A[6,37]) (A[5,37])
post john (`v') (1) (3) (A[1,41]) (A[6,41]) (A[5,41])
post john (`v') (2) (3) (A[1,45]) (A[6,45]) (A[5,45])
post john (`v') (3) (3) (A[1,49]) (A[6,49]) (A[5,49])
post john (`v') (4) (3) (0) (.) (.)

}



postclose john

cap postclose john
postfile john T timing ses beta upper lower using "$table\PreDeath_2", replace

forvalues v=2/4 {
areg death ib4.t2##SES_`v' male i.age married non_dk if low_ses==0 & inrange(t2,0,4) & age0<64, cluster(ydernr) a(ydernr)
mat A=r(table)
mat list A

post john (`v') (0) (0) (A[1,9]) (A[6,9]) (A[5,9])
post john (`v') (1) (0) (A[1,11]) (A[6,11]) (A[5,11])
post john (`v') (2) (0) (A[1,13]) (A[6,13]) (A[5,13])
post john (`v') (3) (0) (A[1,15]) (A[6,15]) (A[5,15])
post john (`v') (4) (0) (0) (.) (.)


areg death ib4.t2##SES_`v' male i.age married non_dk if low_ses==1 & inrange(t2,0,4) & age0<64, cluster(ydernr) a(ydernr)
mat A=r(table)
mat list A

post john (`v') (0) (1) (A[1,9]) (A[6,9]) (A[5,9])
post john (`v') (1) (1) (A[1,11]) (A[6,11]) (A[5,11])
post john (`v') (2) (1) (A[1,13]) (A[6,13]) (A[5,13])
post john (`v') (3) (1) (A[1,15]) (A[6,15]) (A[5,15])
post john (`v') (4) (1) (0) (.) (.)


areg death ib4.t2##SES_`v'##low_ses male i.age married non_dk if inrange(t2,0,4) & age0<64, cluster(ydernr) a(ydernr)
mat A=r(table)
mat list A

post john (`v') (0) (3) (A[1,37]) (A[6,37]) (A[5,37])
post john (`v') (1) (3) (A[1,41]) (A[6,41]) (A[5,41])
post john (`v') (2) (3) (A[1,45]) (A[6,45]) (A[5,45])
post john (`v') (3) (3) (A[1,49]) (A[6,49]) (A[5,49])
post john (`v') (4) (3) (0) (.) (.)

}


postclose john



preserve 
clear 
use "$table\PreDeath_2"
replace timing=timing-4

replace beta=beta*100
replace upper=upper*100
replace lower=lower*100

replace upper=0 if timing==0
replace lower=0 if timing==0



tw 	(line beta timing if ses==0 & T==3, color(navy%60)) ///
	(line beta timing if ses==0 & T==4, color(navy)) ///
	(line beta timing if ses==1 & T==3, color(sienna%60)) ///
	(line beta timing if ses==1 & T==4, color(sienna) ///
	legend(order(1 "High SES Q>=3" 2 "High SES Q>=4" ///
		3 "Low SES Q>=3" 4 "Low SES Q>=4")  region(style(none)) rows(2) pos(6)) ///
	yscale(range(-0.4(0.2)0.4)) ylabel(-0.4(0.2)0.4) ///
	graphregion(color(white) ilcolor(white) lcolor(white)) bgcolor(white) ///
	ytitle("Percentage points") xtitle("Years since clinic closure")	)
graph export "$fig\Death_placebo.png", replace 
graph export "$fig\Death_placebo.pdf", replace 


tw	(rcap lower upper timing if ses==0 & T==4, color(navy%50) lc(*0.5)) ///
	(line beta timing if ses==0 & T==4, color(navy)) ///
	(rcap lower upper timing if ses==1 & T==4, color(sienna%50) lc(*0.5)) ///
	(line beta timing if ses==1 & T==4, color(sienna) ///
	legend(order(2 "High SES" 4 "Low SES") region(style(none)) rows(1) pos(6)) ///
	yscale(range(-0.4(0.2)0.4)) ylabel(-0.4(0.2)0.4) ///
	graphregion(color(white) ilcolor(white) lcolor(white)) bgcolor(white))
graph export "$fig\Death_placebo_Q5.png", replace 
graph export "$fig\Death_placebo_Q5.pdf", replace 

restore


